subroutine rhs(uu,Q)    !UU for Uij,Q for Qrhs
use var
implicit none

    real(8) uu(1+Klx:gridsx+krx,1+Kly:gridsy+kry,5)
    real(8) Q(1:gridsx,1:gridsy,5)
    integer i,j,k
!$omp parallel private(i,j,k)
    !$omp do
    do i=1,gridsx
        do j=1,gridsy
            do k=1,5
                Q(i,j,k)=(flux1(i,j,k)-flux1(i-1,j,k)+&
                         flux2(i,j,k)-flux2(i,j-1,k))/volume(i,j)!(flux1(i,j,k)-flux1(i-1,j,k))/dx+ (flux2(i,j,k)-flux2(i,j-1,k))/dy
            enddo
            Q(i,j,2)=Q(i,j,2)-uu(i,j,1)*gx
            Q(i,j,3)=Q(i,j,3)-uu(i,j,1)*gy
            Q(i,j,4)=Q(i,j,4)-uu(i,j,2)*gx-uu(i,j,3)*gy
            !Q(i,j,5)=(gamma1(i,j,1)-gamma1(i-1,j,1))/dx+ (gamma2(i,j,1)-gamma2(i,j-1,1))/dy
            !write(*,'(I4,I4,2X,100(ES12.4,1X))') i,j,flux1(i,j,3),flux2(i,j,3)
            !print *,i
        enddo
    enddo
    !do i=1,gridsx
    !    do j=1,gridsy
    !        do k=1,5
    !            Q(i,j,k)=(flux1(i,j,k)-flux1(i-1,j,k)+&
    !                     flux2(i,j,k)-flux2(i-1,j,k))/volume(i,j)!(flux1(i,j,k)-flux1(i-1,j,k))/dx+ (flux2(i,j,k)-flux2(i,j-1,k))/dy
    !        enddo
    !        Q(i,j,2)=Q(i,j,2)-uu(i,j,1)*gx
    !        Q(i,j,3)=Q(i,j,3)-uu(i,j,1)*gy
    !        Q(i,j,4)=Q(i,j,4)-uu(i,j,2)*gx-uu(i,j,3)*gy
    !        write(*,'(I4,I4,2X,100(ES12.4,1X))') i,j,flux1(i,j,3),flux2(i,j,3)
    !        print *,i
    !    enddo
    !enddo
    !$omp end do
!$omp end parallel
end subroutine rhs